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Abstract. The LARES satellite, a laser-ranged space experiment to contribute to geophysics observation, 
and to measure the general relativistic Lense-Thirring effect, has been observed to undergo an anomalous 
along-track orbital acceleration of -0.4 pm/s 2 (pm := picometer). This thermal “drag” is not surprising; 
along track thermal drag has previously been observed with the related LAGEOS satellites (-3.4 pm/s 2 ). It 
is hypothesized that the thermal drag is principally due to anisotropic thermal radiation from the satellite’s 
exterior. 

We report the results of numerical computations of the along-track orbital decay of the LARES satellite 
during the first 126 days after launch. The results depend to a significant degree on the visual and IR 
absorbance a and emissivity e of the fused silica Cube Corner Reflectors. We present results for two values 
of am = sir: 0.82, a standard number for “clean” fused silica; and 0.60, a possible value for silica with 
slight surface contamination subjected to the space environment. 

The heating and the resultant along-track acceleration depend on the plane of the orbit, the sun position, 
and in particular on the occurrence of eclipses, all of which are functions of time. Thus we compute the 
thermal drag for specific days. We compare our model to observational data, available for a 120-day period 
starting with the 7th day after launch, which shows the average acceleration of -0.4 pm/s 2 . With our model 
the average along-track thermal drag over this 120-day period for CCR am = tm = 0.82 was computed 
to be -0.59 pm/s 2 . For CCR am = em = 0.60 we compute -0.36 pm/s 2 . 

LARES consists of a solid spherical tungsten sphere, into which the CCRs are set in colatitude circles. 
Our calculation models the satellite as 93 isothermal elements: the tungsten part, and each of the 92 Cube 
Corner Reflectors. The satellite is heated from two sources: sunlight and Earth’s infrared (IR) radiation. 

We work in the fast-spin regime, where CCRs with the same colatitude can be taken to have the same 
temperature. Since all temperature variations (temporal or spatial) are expected to be small, we linearize 
the Stefan-Boltzmann law and, taking advantage of the linearity, we perform a Fourier series analysis. The 
variations are indeed small, validating our Fourier analysis. 

PACS. 00 General physics - 04.80.Cc Experimental tests of gravitational theories 


Contents 


1 Introduction. [5] 

2 Satellite structure . [3] 

3 A few simple estimations . [3] 

4 CCR view factors. [5] 

5 Orbital geometry. [S] 

6 CCR heating and thermal drag. [7] 

7 Computation of the temperatures and the thermal drag. [9] 

8 Conclusion . ED 

9 Acknowledgement . [TO] 

A Radiative heat exchange. |T0] 

B Numerical values of the constants. m 

C Analytical computation 1: the temperature of a stationary spinning sphere bathed in sunlight. [19] 

D Analytical computation 2: Going beyond the linearized boundary condition. [23] 


a Email: phn229@physics.utexas.edu 
b Email: matzner2@physics.utexas.edu 

















2 


Phuc H. Nguyen, Richard Matzner: Modelling LARES temperature distribution and thermal drag 


1 Introduction 


LARES is a space mission sponsored by the Italian and European space agencies, which consists of a passive satellite 
tracked by the global network of laser ranging stations of the International Laser Ranging Service. The satellite is a 
ball of tungsten 18.2 cm in radius and weighing about 400 kg, with 92 fused silica ( Suprasil ) Cube Corner Reflectors 
(CCRs) covering the surface. The orbit is nearly circular, with a perigee of 1450 km and an orbital inclination of 
69.5 degree. It is similar in structure (but smaller and denser) to the related satellite LAGEOS. LARES consists of 
a solid tungsten core, whereas the LAGEOS satellites have a cylindrical brass core surrounded by an aluminum shell 
consisting of two joined hemispheres, into which the CCRs are set. 

One of LARES’s scientific targets is to verify the Lense-Thirring effect [HIIDEj (also known as frame-dragging), a 
prediction of the theory of General Relativity, using the technique of laser ranging with an accuracy of about 1 % [H 
rnrnm ■ Frame-dragging is the phenomenon by which a spinning mass drags the local inertial frame in the direction 
of the rotation. This induces a precession of the plane of an orbiting test body ;9j. The ascending node of an orbit 
around the Earth drifts eastward. In the weak field slow rotation limit appropriate to orbits around the Earth, this 
precession is given by given by: 


fl 


Lense-Thirring 


2 J 

o 3 (l - e 2 ) 3 / 2 


( 1 ) 


For an object in low orbit around Earth, in particular, the formula above yields a precession of the order of a tens or 
hundreds of milliarcsecond per year, corresponding to meters per year motion of the node. Studies using Earth gravity 
field from the GRACE mission together with the laser ranged satellites LAGEOS and LAGEOS 2 have observed the 
node dragging to approximately 10 % .10]. Due to the smallness of the Lense-Thirring effect, experimental detection 
is a real challenge and it is crucial to be able to model all the experimental sources of errors as accurately as possible 
□ 

There are many contributions to the nodal precession rate, both gravitational and non-gravitational. Many of these 
perturbations are orders of magnitudes larger than the Lense-Thirring effect and not all of them are well-modelled. 
Precession due to Newtonian effects arises from the fact that Earth is not a perfect sphere. Non-gravitational effects 
include thermal and residual-gas effects. 

In the rest of this paper, we will focus on the measure of one of the most important non-gravitational sources of error: 
thermal drag. This is the net deceleration that the satellite undergoes due to a nonuniform temperature distribution 
(and therefore a nonuniform thermal radiation at the surface). In the case of LAGEOS, this net force explains that 
satellite’s anomalous along-track acceleration of about -3.4 cm/s 1 2 . For LARES’s first 126 days in orbit, its average 
drag over the last 120-day period was measured as -0.4 pm/s 2 [12j. We will model the temperature of the LARES 
satellite with the following two assumptions (justified below): (1) the temperature variation inside the metal part of 
the satellite or within a single CCR is ignorable compared to the temperature variation between the metal and any 
CCR or between different CCRs, and (2) the satellite spin is sufficiently fast that CCRs of the same colatitude (with 
respect to the spin axis) have the same temperature. 

The rest of the paper is organized as follows. In Section [2j we describe the structure of the satellite in detail. In 
Section [3j we estimate in a simple way the temperature variation within the tungsten, and within each CCR, thus 
justifying the assumption that temperature variations inside each CCR and the metal are ignorable. In Section [4j we 
compute the so-called view factors of the CCR cavity, which are necessary ingredients to handle the radiative heat 
transfer inside the cavity. In Section [5] we describe the geometry of the orbital plane, the spin orientation, and the 
solar eclipse. In Section [6j we give an overview of the satellite’s thermal response and thermal drag. In Section [7J we 
present the results of our calculation. In Section [8j we summarize the main findings of our study. In Appendix [Aj 
we review the main elements needed to handle radiative heat exchange problems, in particular the formalism of view 
factors and the radiosity method. In Appendix[B] we list the numerical values of all the constants used throughout the 
paper. Finally, in the two last appendices (Appendix [C] and |D[) , we present two analytical computations to find the 
temperature inside a metal sphere bathed in sunlight /"Appendix 0 focuses on the case of a spinning sphere without 
orbital motion. While these computations are not directly related to the main points of the paper, they employ several 
of the tools also used in the paper (especially the Fourier series technique), and can be of interest to future researchers 


1 The frame dragging leads also to a related effect, precession of the direction of a gyroscope: the direction of the gyroscope 
which initially points to a specific celestial coordinate (eg toward a particular QSO), deviates from this direction at the rate 


GI 

W 9V ro ~ Z^/3 




( 2 ) 


where uj gyro is the gyroscope precession, / is the moment of inertia of Earth, r is the position of the gyroscope, and S is Earth’s 
angular velocity. At the 642 km altitude of the GP-B experiment m, equation pi gives a precession rate of —39.2 milliarc¬ 
second per year, where the minus sign indicates that the precession is Eastward ( H] Fig. 1). GPB observed a frame-dragging 
of —37.2 ± 7.2 milliarc second per year, about a 19% quoted error. 
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Fig. 1 . The shape of the CCR. 


in the field of radiative heat exchange. Finally, Appendix [D] considers a stationary, nonspinning sphere, and discusses 
a possible treatment to go beyond the linearized boundary conditions. 


2 Satellite structure 

In this section, we describe the arrangement of the CCRs on the satellite, and the shape of each CCRs [TiTin~T| . There 
are 92 CCRs on LARES, and they are arranged into 10 rows. We give below the number of CCRs per row (Nj where 
I labels the row) and the colatitude 9j of the row. The CCRs belonging to a given row are equally spaced. 


Row 

ni 

0i 

Row 

ni 

9i 

I 

1 

0 deg 

-V 

16 

100 deg 

II 

5 

20 deg 

-IV 

14 

120 deg 

III 

10 

40 deg 

-III 

10 

140 deg 

IV 

14 

60 deg 

-II 

5 

160 deg 

V 

16 

80 deg 

-I 

1 

180 deg 


The shape of each CCR is depicted in fig. |T| The CCR is a right triangular pyramid intersected by a cylinder of 
radius R. The CCR fits inside a cylindrical tungsten cavity with a conical bottom. We will denote by d (« 5 mm) the 
distance between the tip of the CCR and the bottom of the conical floor of the cavity. For the sake of visualization, 
we have labelled a few points on fig. [T] and given the coordinates of those points in the table below. The origin of the 
coordinate system was chosen to be at the tip of the CCR. 
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3 A few simple estimations 

In this section, we justify the assumption of uniform temperature inside the metal and each CCR by estimating the 
temperature variation within the metal and within a CCR. First suppose the satellite to be a simple tungsten ball 
in radiative exchange with empty space, and heated by sunlight (in this section we will ignore Earth IR radiation). 
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To compute the power absorbed by the ball, notice that the heat absorbed by each area element dA on the sphere is 
proportional to the cosine of the angle between the normal to dA, i.e. the cosine of the colatitude of dA. Integrating 
over the illuminated hemisphere then yields: 

p2n /♦ 7 t /2 

Q = aw,vis&Rt a t / / cos 9 sin 9ddd(f> (3) 

Jo Jo 

where R sa t is the radius of the LARES satellite, aw,vis is the absorptivity of tungsten in the visible and d> is the solar 
constant. The values of all constants appearing in this paper are tabulated in Appendix [B] Evaluating the integral 
above, we find: 

Q = a WtVis <PnR 2 sat (4) 

i.e. the heat absorbed is proportional to the cross-section nR 2 at . 

By Fourier’s law of heat conduction, Q is roughly the product of the temperature gradient, the cross-section area 
through which heat is conducted, and the thermal conductivity Kw of tungsten: 


AT 2 

( 5 ) 

Q ~ sat 

P^sat 

The two equations above yield 


^ rj-, ^ ^W,vis^Rsat ^ 

(6) 

Kw 

On the other hand, conservation of energy reads: 


ew,iR<rTw4nR% at = aw,vis^TrR 2 * at , 

( 7 ) 


where ew,iR is the emissivity of tungsten in the IR, a is the Stefan-Boltzmann constant, and Tw is the metal’s 
temperature. From this equation, we find: 


T w =( °" w ’ vis$ \ ' « 443.6K (8) 

\4 zw,irg ) 

Thus, the temperature variation AT is only about 0.2 percent of the average temperature. We can therefore take the 
metal to be at a uniform temperature to an excellent degree of accuracy. 

Next, we turn to the estimation of the temperature variation inside a CCR. For simplicity, we can model the CCR as 
a cone with half-angle j inside a cylindrical cavity of radius and height R, and such that the base of the cone is flush 
with the satellite’s surface. The CCR is in radiative heat exchange with the metal wall of the cavity, which is assumed 
to be at constant temperature TW as previously found. We will assume no conduction takes place from the metal to 

theCCRsEl 

An exercise in radiative heat exchange (see Appendix [A| shows that the net power transferred from the metal to the 
CCR is given by: 

Pnet = e effAgia(T w — T gl ) (9) 

where A g i = V2 (ttR 2 * ) is the surface area of the glass wall, T g i is the average CCR temperature, and e e // is an effective 
emissivity which is solely determined by e g i,iR , cwjr as we ll as the geometry of the cavity. To compute e e ff, we need 
the so-called view factors between isothermal elements lining the cavity F^j , which represent the fraction of radiation 
leaving isothermal element i which strikes element j. In our case, we only have 2 isothermal elements: the glass and 
the metal. Moreover, since the CCR is a convex surface, two of the view factors are trivial: 

Fgi.gi = 0 ( 10 ) 

and 

F g i,w = 1 (11) 

The other two view factors follow from the surface area of the glass (A g i) and the metal {Aw) lining the cavity: 

A gl = V2i tR 2 (12) 

2 In reality, the CCRs are not in contact with the metal, but are held in place by Teflon spacers, which are covered at the 

surface by a tungsten retainer ring. Teflon is a very good thermal insulator, so we assume the CCRs have only radiative contact 

with the metal. Moreover, we will ignore the effect of the retainer rings; the temperature of a copy of LAGEOS 2’s retainer 

rings was measured in a 2006 laboratory test. The retainer rings were found to have closely the same temperature as the body 

of the satellite m- Hence the retainer rings should not have any noticeable effect on the thermal drag. 
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A\y = 3nR 2 


Using the identities (74) and (75) for view factors, we find: 


F W ,„ = f 


and 


Fw,w = 1 — 


V2 


Finally, the effective emissivity is given in terms of the view factors by: 


Specifically, in our case: 


1 1 - e w ,iR „ 

£eff = -1- *W,gl 

£gl,IR ew,IR 

1 , 1 — tw. l Ft 'JF' 1 

e eff ~ ' 


TgUR e w,iR. 3 

with tgi.iii the emissivity of glass in the IR. This net heat is then radiated into space by the base of the cone, which 
faces into space: 

Pnet = tglJRO’TgiirR 2 


We can now solve for the average glass temperature [^] 


Tgl = ( 1 


£gl,IR 


- 1/4 


T w « 291.8 K 


V^Ceff, 

Using Fourier’s Law, we can approximate the heat flow in the glass just below the surface to: 

„ AT n2 

Pnet — Kgi TtR 

where n g i is the thermal conductivity of glass. The temperature variation inside the CCR is then: 

teffAgMTL-T*) 

AT = 11 y w -« 3.84K 

KginR 

This represents 1.3 percent of the average CCR temperature. This result is consistent with experimental (ground 
based) measurements of the CCR temperature distribution [ 15 ] . 


(13) 

(14) 

(15) 

(16) 

(17) 
/hich 

(18) 

(19) 

( 20 ) 

( 21 ) 


4 CCR view factors 

We will assume that only radiative heat transfer takes place between the CCR and the tungsten, and no heat is 
conducted. In order to use the radiosity method, we will need to write down the view factors and the effective 
emissivity for the geometry depicted in fig. [T] First, the total metal surface area lining the cavity is: 

A w = 2ti -R{y/2R - 2d) + ttR^/R 2 + 9 d 2 (22) 

The total glass surface area lining the cavity is: 

A gl = (V37T + 2v^tt - 3 V6)R 2 (23) 

Notice that, like in Section[3] the CCR is still a convex surface, implying that two of the four view factors still vanish: 

Fgi,w — 1 (24) 

Fgi,gi = o (25) 


The other two view factors can be found using the areas Aw and A g i, and the relations (74) and (75): 

(3\/6 — \/3n)R 2 — 47ri?.d + ttR\/ R 2 + 9 d 2 


Fw,w = 


Fw,gi = 


2nR(V2R - 2 d) + nRVR 2 + 9 d 2 

(\/37r + 2 y/2-n - 3 V6)R 2 

2itR(V 2R - 2d) + nRVR 2 + 9 d 2 


(26) 

(27) 


3 


We are assuming that the CCR is not illuminated by the Sun in this simple example. 
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5 Orbital geometry 

In this section, we describe the geometry of the orbital motion and the spin axis (see fig. [2] for an illustration of the 
orbital geometry on the day of launch). The orbit is nearly polar, highly circular with inclination i = 69.5 degrees, 
which we take as 70 degrees in our model. Due to the Earth’s quadrupole moment, the orbital plane precesses westward 
at a rate of 1.7 degree/day. This value, and all values for the orbital parameters for the first 126 days are taken from 
M- The satellite was successfully launched in orbit on February 13, 2012 at 10:00 UTC from the European Space 
Agency’s spacesport in Kourou, French Guyana m- To describe the orbit, we will work in celestial coordinate system 
2 The unit vector rg U n pointing from Earth to Sun k days after launch is: 

/ cos (27r(fc - 37)/365) \ 

rsun(k) = I cos/3sin (27r(fc — 37)/365) I (28) 

ysin/?sin ( 27 t (k — 37)/365) J 

where /3 = 23.2 degree is the obliquity of the Sun’s orbital plane (i.e. the angle this plane makes with the Earth 
equatorial plane), and 37 is the number of days from the date of launch of LARES to the vernal equinox. We will 
assume that for any given day, the position of the Sun in the sky and 17 stay fixed (i.e. they increment discretely from 
one day to the next). Also, the unit vector r sat from the center of Earth to the center of the satellite is: 

( — sin 17 cos i sin oj 0 t + cos 17 cos ui 0 t\ 

cos 17 cos i sin ui 0 t + sin Q cos u> 0 t (29) 

sin i sin u> 0 t J 

where wo is the satellite’s orbital frequency, and 

17(fc) = (220- 1.7fc) I ^ (30) 

is the longitude of the ascending node on day fc[^] The time t will be taken to range in t £ [0, T], i.e. one orbital period, 
with the satellite at the longitude of ascending node at t = 0. 

The spin orientation and spin frequency have been measured and the results are reported in [18| . The spin orientation 
is at RA = 12 ft, 22 m 48 s ( RMS = 49 m ) and Dec = —70.4 degrees ( RMS = 5.2 degrees). For simplicity, we will take 
the spin direction to be (in celestial coordinates): 


S = 



(31) 


where i = 70 degrees is the orbital inclination of LARES. The above choice for S is within the margin of error for the 
RA and Dec. Next, the unit vector fccR pointing from LARES’s center to a CCR at colatitude dj is: 

( — sin i sin 0/ cos uit — cos i cos 0 A 

sin 8j smut (32) 

cos i sin 0 / cos ujt — sin i cos 0 / I 

where w is the spin frequency. Notice that the time dependence in the vector above is entirely due to the spinning of 
the satellite. The spin frequency actually decreases over time due to interaction with the Earth’s magnetic field, as 

am): 

uj(k) = (0.546rad/s) exp -0.00322509A: (33) 

For k = 0, the spin angular frequency is about 600 times the orbital frequency. For k « 100, the spin frequency has 
decreased to about 400 times the orbital frequency. Thus, we can see that during the first 126 days, we remain safely 
inside the fast-spin regime. 

4 i.e. Cartesian coordinate system, with the Earth’s equatorial plane as the x-y plane, and the x-axis pointing toward the 
vernal equinox. 

5 For a geocentric orbit such as LARES, the longitude of the ascending node (also known as the right ascension of the 
ascending node or RAAN) is the angle, measured eastward, from the vernal equinox to the ascending node. The ascending node 
is defined to be the intersection of LARES’s orbit with the Earth’s equatorial plane, where LARES is heading north. 
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Fig. 2. The situation at t = 0 and k = 0, as viewed directly above the North Pole of Earth. The sphere centered at the origin 
represents the Earth. The coordinate axes form the celestical coordinate system. The orange sphere is the Sun. The dashed gray 
circle is the Sun orbit through 1 year. The gray sphere is LARES, with the vector representing the spin orientation. The blue 
ellipse is the projection of LARESs trajectory for k = 0 (i.e. day zero, launch day) on the x-y plane. The radial distances from 
the Earth’s center are not drawn to scale. However all angular information in the sky is drawn accurately. 


6 CCR heating and thermal drag 

The physical concepts that lead to along-track thermal drag for LAGEOS and LARES style satellite (passive metal 
satellites with CCR retroreflectors) were developed over a number of years by |T 9 II 20 y 2 lll 22 y 231 l 241 l 25 [[ 2 i 3 lj and incor¬ 
porated into Slabinski’s definitive study of the along track thermal drag on LAGEOS m- Our analysis is based on 
these references, though with technical differences; we do Fourier analyses while Slabinski did direct time integration, 
for instance. In this Section we discuss the two concepts most crucial to understand the LARES along-track thermal 
drag. Our approach will be to linearize the thermal equations and develop our study in terms of Fourier analysis. This 
means that the heating due to Earth IR and that due to solar heating can be treated, at least conceptually, separately. 
In our actual computations we sum all contributions. 

As we have noted, the tungsten core of LARES has very high conductivity, and small temperature differences (which 
we ignore) across its surface, regardless of the heating source. However the CCRs are glass which is a good heat 
insulator, has a substantial heat capacity, and is an efficient infrared radiator. These properties mean that the CCRs 
have substantial thermal inertia. Thus the CCRs heat gradually when exposed to a heating source, and cool gradually 
when the source is removed. 


6.1 Earth Infrared CCR heating and thermal drag 

The thermal drag from Earth-IR heating is due to the Earth Yarkovski effect [22] (see fig. 3] for an illustration of this 
effect). The simplest example has the spin axis in plane of orbit. Maximum hemispheric IR absorption occurs when 
one pole points toward Earth IR source. Because the CCR -glass- is opaque to IR, with low thermal conductivity, 
CCRs near the heated pole warm slowly, so the along axis thrust builds up after the satellite has moved from the the 
position of maximum hemispheric IR absorption. (In the fast spin case we consider, the temperature distribution is 
axisymmetric around the spin axis.) The along-track thermal drag (it is always a drag) maximizes after some phase 
angle in the orbit. The maximum effect occurs twice per orbit. 

If the spin axis is out of orbital plane there is less differential heating, and the need to project resulting (smaller) 
thrust along orbit reduces the along-orbit thermal drag further. But it is always a drag. 
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Fig. 3. Illustration of the Yarkovsky effect due to Earth IR radiation. The satellite is depicted at 8 different positions along 
its orbit, and its spin axis is the arrow in magenta at the bottom. The direction of motion of the satellite is clockwise. The 
hotter hemisphere is depicted in orange, and the cooler hemisphere is depicted in blue. A sharper color contrast between the 
two hemisphere indicates a greater thermal anisotropy. When there is no anisotropy, the satellite is depicted in gray. The length 
of the black arrow indicates the magnitude of the axial force. For simplicity, we depict the situation where the thermal lag phase 
angle is 7r/4. 


6.2 The eclipse 

For certain configurations of orbital plane and sun position, the satellite will undergo (solar) eclipses as it passes 
behind the Earth, i.e. into the shadow of the Earth, and is temporarily shielded from direct solar radiation. 

When there are no eclipses, in the fast spin regime solar radiation has no effect on the thermal drag averaged over 
an orbit. This is because the spin is parallel propagated, so always points to the same celestial coordinates, and for 
one orbit, the solar celestial position is very closely constant also. Thus the same hemisphere of the satellite is heated 
throughout the entire orbit. The consequent spin-axis thrust is constant in size, and constant in direction (along 
the fixed spin axis) throughout the orbit. However, when dotted into the orbital velocity (to obtain the along-track 
acceleration), it is alternately a positive acceleration, and a drag on the opposite side of the orbit, so that the along 
track thermal drag arising from this source averages out over one orbit. 

When there is an eclipse, however, the heating is removed for part of the orbit. If we incorrectly assume that the 
thermal thrust immediately disappears during the eclipse, then the along track acceleration during the eclipse is no 
longer available to cancel the contribution from the other side of the orbit. This leads to a nonzero net along-track 
acceleration. This net acceleration can be either positive, or negative. In his LAGEOS studies, Slabinski [27] found 
that the size of the effect could be as large as the Earth-IR induced thermal drag, and because of the possibility of 
either sign depending on the spin/orbital plane/sun geometry, could swing the thermal drag to twice the noneclipse 
thermal drag, or zero, consistent with LAGEOS observations. 

In both LARES and LAGEOS, the thermal inertia of the system means that the thermal thrust does not disappear 
immediately upon entering the eclipse, but decays slowly, and recovers gradually when the satellite emerges from 
the eclipse. But the principle is as explained in the previous paragraph. Also note that because of the lower orbit of 
LARES, it spends a much larger fraction of its time in eclipse than does LAGEOS. 

We will need to compute the time of eclipse entrance t en t ra nce and eclipse exit t ex u . To do this, we compute the 
perpendicular distance d(t, k) from the satellite to the Earth-Sun axis: 

d(t, k) = a\v sa t (t, k ) x r Su n{k)\ (34) 

Then t entra nce and t ex u are solutions to the equation (see footnote [6] below): 

d(t, k) = Re 


(35) 
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However, in general there will be 4 solutions. To see this, one can visualize a cylinder with the Earth-Sun axis for axis 
of symmetry, and with the radius of Earth for radius. The part of the cylinder “behind” the Earth is in shadow. Thus 
if, on a particular day, LARES experiences an eclipse, then its orbit will intersect the cylinder in 4 points, 2 of which 
are t en t ra nce and t ex u . To pick out the correct two solutions, we consider the dot product rsun ■ f sa t- When this dot 
product is negative, the vector fsun makes an angle larger than ^ with the vector f sa t, and there is eclipse. 


7 Computation of the temperatures and the thermal drag 

In this section we will compute the temperature of each isothermal element, and subsequently the thermal drag. 


7.1 Surface condition of the CCRs 

Appendix |B| lists parameters for our model, including the absorptivity/emissivity of the fused silica CCRs both in the 
visible andin the IR. They are quoted both for clean glass - measured in lab, and for “on-orbit”, or “dirty” glass. 
Absorptivity and emissivity are surface properties, and are modified by surface contamination. The on-orbit visual 
and IR values are important because they determine the heat balance in the CCRs. 

To begin with the values in the visual, clean Suprasil glass has a very small absorptivity in the visual (a few parts per 
million [25]), and is quite dark in the IR (em = am = 0.82). As Slabinski [27] pointed out, an aside to discuss Optical 
Solar Reflectors (OSRs) is useful at this point. OSRs are silvered glass mirrors which are placed (silver side in) on 
spacecraft exteriors to control temperature in satellites exposed to sunlight. The silver reflects most of the incident 
visual energy, but the glass, which has a high emissivity in the IR, efficiently radiates heat away. 

In the 1970s it was observed that OSRs became less effective over time. At least two experiments [ 1130 ] were flown 
with a carefully designed and modeled configuration of heater, thermocouples, and OSRs, to measure the visual 
absorptivity of the OSRs. It was found that the visual absorptivity increased from about 10 % after launch to over 
15 % over a period of approximately a year. In his exhautive paper on thermal effects in the LAGEOS satellite, 
Slabinski m took a value of a V i S = 0.15. He quotes conversations with aerospace designers that “surfaces in the space 
environment tend toward gray. Initially white surfaces turn to gray; initially black surfaces turn to gray” (Slabinski, 
private communication 2015). We follow Slabinski by taking the on-orbit = 0.15 for the CCRs. 

However, Slabinski uses a clean glass IR emissivity: 0.82. We show below that the somewhat smaller value of 0.60 
matches the LARES observations more closely. Physically, any modification to the glass surface should change the IR 
properties as well as the visible properties, and moving the IR emissivity 0.82 —> 0.60 moves in the direction of graying 
the dark IR surface of the clean glass. Thus we regard an on-orbit value of 0.60 a plausible value of IR emissivity for 
the CCRs. Reducing em = am reduces the thermal response of the satellite to warming by the Earth, which reduces 
the thermal drag. We present results for both em = am = 0.82, and for em = am = 0.60, below. 


7.2 The radiation intensities 

First we will need the radiation intensity incident on each CCR. The contribution from sunlight is quite straightforward. 
First we compute the angle between the normal direction to the CCR and the sunlight direction: 

cos (7 vis (t, 9 ] r , k)) = rccR(t, 0i) ■ fsun(k) (36) 

Notice that we are assuming the Sun to be a point source, so that all the sunlight comes from the direction along the 
Earth-Sun axis, since the Sun subtends a very small solid angle in the sky, as seen from LARES. The visible intensity 
incident on the CCR is then: 


Ivis{t,9i,k) = &cos (y via (t, 9/, k))&(cos (y vi 3 (t, 9j, k))) (37) 

where the unit step function O ensures that, when the CCR faces away from sunlight, the intensity vanishes (at any 
moment only half the satellite is illuminated by sunlight). Moreover the above is only valid outside the solar eclipse; 
the visible intensity also vanishes when the satellite is in the eclipse. 

Next, consider the Earth IR radiation. Here we hit a substantial complication owing to the fact that the Earth, unlike 
the Sun, cannot be modelled as a point source. We have to take into account the finite-size effects of the Earth. 
Fortunately there is a remarkable simplification if we assume the Earth to be a Lambertian emitter: as viewed from a 
unit surface area of the satellite, the Earth’s intensity per unit solid angle is a constant, given by Nm = 71W• m -2 • sr -1 
([27]). More explicitly, this means that for each unit of solid angle of the source (i.e. the IR Earth) as seen from the 
satellite, a unit area of the receiver’s surface receives 71 Watts if oriented perpendicularly to the radiation (if the 
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Fig. 4. The setup used to compute the Earth’s IR radiation on the satellite Iir(9ccr), with the larger sphere representing 
Earth and the smaller sphere representing LARES. Here we have depicted a situation with 9ccr > 0 (i.e. the CCR faces toward 
Earth). The portion of the celestial sphere subtended by the Earth from LARES, over which we integrate to obtain Iir, is the 
spherical cap, which has an angular radius 54.55 degrees. The separation of LARES and the Earth is not to scale in the figure. 


orientation is not perpendicular, the power will be reduced by a Lambertian cosine factor, as usual). Thus, the total 
IR intensity is obtained by evaluating an integral of the solid angle subtended by Earth as viewed from LARES, 
weighted by the Lambertian cosine factor, then multiplying the result by Nir. 

To evaluate the integral described above, we set up spherical coordinate system (r, 0,0) centered on the satellite such 
that the z-axis is along the satellite-Earth axis (see fig. [4]) . Then the Earth subtends the disk 9 < a e = 54.55 degrees, 
where 


a e = arcsin 



(38) 


where Re is the radius of Earth’s IR-emitting atmosphere]^] Moreover, notice that there is no special direction for the 
spin axis in these coordinates. Thus we can without loss of generality rotate coordinates around the z-axis to assume 
that the CCR momentarily lies in the x-z plane. Then the unit vector fccR from the satellite’s center to the CCR 
can be parametrized by one angle Occr (see fig- 13 : 


rcCR = (cos dccRi 0, sin 9 C cr) (39) 

with 9ccr ranging from — ^ to f. Then we can write down the cosine of the angle between fccR and the vector from 
the CCR to some area element on the Earth disk with coordinate (0, 0): 

COS (t ir(9ccr, 0, 0)) = COS 9ccr sin 0 cos 0 + sin 9 C cr cos 0 (40) 

The IR intensity on the CCR is now only a function of 9ccr- It can be piecewise defined as follows. If Occr < —cn e , 
then by inspection of fig. [4] we can easily seen that the CCR does not intercept any IR radiation. In this case 
Iir(9ccr) = 0. If — a e < 9ccr < 0, the CCR sees less than half the Earth disk, and 

poc e rF(0ccR,0 ) 

Iir{9ccr) =2Nir / sin 0 d 0 / cos {"/ir(9ccr,0, 0 ))d 0 (41) 

J\0ccr\ Jo 

6 Technically, for the purpose of the eclipse, we should use the actual Earth radius (6378 km) rather than the radius of the 
IR-emitting atmosphere, which is slightly larger (6407 km). However the difference between the two radii introduces a negligible 


error. 
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/|R 



Fig. 5. Plot of the function Iir(6ccr) in blue. The two vertical red lines correspond to 9ccr = ±a e - For comparison, we have 
also included the intensity profile if the Earth was approximated as a point source (in dashed gray). 


with F(0ccr,9) obtained by setting equation (401 to zero and solving for <j>\ 


F(0ccR,e) = arccos 

Next, if 0 < 9ccr < «e) then the CCR sees more than half of the Earth disk, and 


(42) 


Iir{0ccr) = 2 N ir 
+ 2 


rF(6ccR,Q) 


sin Odd 


\Sccr\ Jo 

\Qccr\ 

sin Odd / cos'ymiOccR, 0,(j>)d</) 


cos (7 ir(6 C cr, 9, </>))# 


(43) 


Finally, when a e < dccR, then the CCR views the whole Earth disk, and the integral can be evaluated analytically: 


Iir{0ccr ) = kNi R sin 9 C c r sin 2 a e ( 44 ) 

The function Iir(6ccr) is plotted in fig. [ 5 ] Finally, Qccr itself is a function of the time t elapsed during one orbit, 
the number of days k since launch, and the colatitude 6j of the CCR: 

0CCR(t, k, Qi) = - arcsin {f C CR(t, 0i) ■ f sat (t, k)) (45) 

We close this subsection by a few remarks on the error that would be introduced if we were modelling the Earth as 
a point source. One can multiply Nir by the total solid angle Aft subtended by Earth as viewed from LARES to get 
an effective power per unit receiver area Wir (of the dimension W • m -2 ). The solid angle AQ is easily computed: 


r&e. 

Afl = 2tt / sin 6d6 = 2 tt (1 — cos a e ) 

J 0 


(46) 


We then find Wir = 188W • m -2 . We must multiply this flux by the cosine of the angle between the CCR normal 
and the direction of the radiation, i.e. cos (| — Occr ) to obtain the power absorbed by a surface element. In fig. ol 
we also plot this resulting intensity profile for comparison to the correct treatment. As can be seen from the plot, 
the finite-size effect of the Earth spreads out the intensity more evenly across the satellite surface, thus significantly 
reducing the thermal drag. 
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7.3 The equations of conservation and the temperatures 


In this subsection, we write down the equations of conservation of energy for each CCR and the metal, Fourier 
transform those equations and solve for their temperatures. For a CCR of colatitude 0/, conservation of energy reads: 

e c „A al „(TUt)-T}{t)) + + (47) 

where Tj is the CCR temperature, which only depends on the row I since we are in the fast-spin regime. The terms 
on the left-hand side represent the following contributions: the first term is the net heat tranferred to the CCR from 
the cavity, the second term is the heat absorbed at the exterior surface from Earth IR radiation, the third term is the 
heat absorbed at the exterior surface from sunlight, and the fourth term is the heat lost into space by radiation. On 
the right hand side is the rate of change of the temperature with time. 

For the metal, conservation of energy reads: 


Pw,visit) + Pwjr — A g ie e ffcr ( 92T(^(f) — Y^ mTf{t) ] — ewjRAw,vac&Tw{t) = mwcw 


dT w (t) 

dt 


(48) 


The first term on the left is the power absorbed by the metal from sunlight. Its time-dependence is solely due to the 
eclipse. When the satellite is outside the eclipse, it is approximately constant and approximately given by: 


Pw,vis ^W, vis^ ^W^vis^^R T 


^(1 Otgl, v is) a gl,z 


nR 2 $ E nj cos 9 1 


(49) 


The power absorbed by the metal from sunlight is in fact a funtion of the orientation of the satellite, and hence of 
time, but we assume that it is approximately constant, and evaluate for a specific configuration, namely when one pole 
of the satellite is directly pointed at the Sun. Then the first term above is the power absorbed as if the satellite were 
a simple metal sphere. The second term above is the power incident on the North Pole CCR which, instead of being 
absorbed, is reflected back. The last term above accounts for the remaining CCRs (other than the North Pole). We 
assume that a fraction 4(1 — a g i jV i S ) of sunlight incident on these CCRs make it into the cavity and ends up absorbed 
by the metal (see mi_ 


The second term in (48) represents the power absorbed by the metal from the IR. It is approximately constant in 
time: 

/* ~ —|— <~V c 

- ,2 a ^d 9 — ewRFnrR 2 ^' ~ 61' 


Pwjr = cwjRPsat 


sin 


/o 


0Iir( 


cwjrttR 2 Y^ n idiR^- 


(50) 


where the first term on the right-hand side represents the power absorbed as if the satellite were a metal sphere, 
and the remaining terms account for the CCRs. Coming back to equation (48), the third term on the left-hand side 


represents the power transferred to the CCRs from the metal, and the 4th term on the left-hand side is the power 
radiated into space. 

Since we expect the time variation in the temperature to be small, we linearize the conservation equations by letting 


Tw{t) = T Wfi + AT w {t) (51) 

Tj(i) = T l o + AT T (t) (52) 

with AT\y << Tw, o and ATj « Tj.q. Then: 

Tw(t) « Tw,o + 4 T^ 0 AT w (t) (53) 

!/(*)« I/, 0 + 4lf, 0 ^2>(t) (54) 

Moreover, thanks to the periodicity of the motion, we can expand the time-dependent pieces of the temperatures into 

Fourier series 0 _ 

AT w (i) = Y, T w,ne iui ° nt (55) 

ATiit) = Y Ti,ne iuiont (56) 

n(zZ,r 17^0 


' Notice that we carry out a Fourier sum, not a Fourier integral. To do this we are implicitly assuming that the spin frequency 
LOs is an integer multiple of the orbital frequency lo 0 , so that all time-dependent quantities are exactly periodic. In the rapid-spin 
regime, this does not cause any difficulty (in particular, no resonances). 
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Decomposing the conservation equations into modes, we find, for the time-independent mode: 

e e ffA g icr(Tw t o ~ lf,o) + e gi,iR 7r R 2 IiR,o(0i) + e g i,visirR 2 Ivis,o(9i) — £ g ijR&irR 2 Tf fi = 0 (57) 

P\V,vis, 0 + Pw,ir — A g ie e /M92^, 0 -E njl?, 0 ) — e w ,iRA w ,vac<jT^ 0 = 0 (58) 

i 

where Pw,vis, o> Iir ,o and I v is,o are the time-averaged Pw,vis, Iir. and I V i S , respectively. All of these three quantities 
can be computed by numerical integration. For the n = 1 mode, we have: 

4 e e ffA g ia(Tyy 0 Tw,i - T? gT/p) + e g ijRnR 2 I I R^(9i) + a g i, vis TrR 2 I vis ,i(9j) - 4:e g ij R airR 2 T? 0 T I}1 = mccRC g iiu 0 T I} i 

(59) 

P\V,vis,l -4A sJ e e// a(92T^ 0 T Wil -^ n iTl 0 T IA ) — i iewjRAw,vacvTwfi r Pw 1 i = ’mw c wi’U 0 T w ,i (60) 

/ 

For the n = 2 mode, we have: 

^ e effA g ia(Tw 0 Tw,2 —TjqTi^) + £ g l,IR^R 2 IlR,2^l) + a gl,vis'KR 2 I v is,2(0l) ~ & g l,IRC r 'KR 2 Tf 0 T/ i2 = 2mcCR c gl^ u} oTl ,2 

(61) 

Pw,vis,2 -4A si e e// a(92T^ 0 T w , 2 -^ niTifiTi^) — 4.e W jRA w ,vacvTw,o T w,2 = 2mwcwiu 0 Tw,2 (62) 

i 

By solving the equations above, we obtain the coefficients Tw, o, Tw, i> Tw ,2 and Tj^, T/p, Tj ) 2 . The temperatures are 
then: 

Tw(t) = Tw ,0 + 2Re(2V ) i) coswof — 2Im(IV i i) sinwot + 2Re(7V, 2 ) cos2wof — 2Im(2V i2 ) sin2wot (63) 

Ti(t) = Ti,o + 2Re(T/,i) cosujot — 2Im(T/ i i) sinwot + 2Re(T/ ;2 ) cos2wof — 2Im(T/ i2 ) sin2wot (64) 


7.4 The thermal drag 

We assume the surface of the satellite to be Lambertian: the angular dependence of the intensity emitted by any 
surface element 1(9) is proportional to the cosine of the angle between the normal vector to the surface element and 
the observer’s line of sight: 

1(9) = /o cos 9 (65) 

where Iq is the intensity emitted in the normal direction. Integrating over all angles yields the total radiated intensity 
by the surface element, given by the blackbody radiation formula: 


eaT 4 = 


Iq cos 9dfi = ttIq 


( 66 ) 


where e is the emissivity of the surface element, a is the Stefan-Boltzmann constant, and T is the temperature of the 
surface element. The total momentum emitted in a unit time by the surface element into an angular ring dfi is: 


dP(9) = ^dn 
c 

The force on the surface element is obtained by integrating over the hemisphere: 

F = — [ cos 9dP(9) = — --- 

./ 3 c 


(67) 


( 68 ) 


where, by symmetry, it is enough to integrate over the normal component of dP. We emphasize that this expression, 
in particular the temperature T appearing in equation |68[ refers to one infinitesimal surface ring of the satellite; the 
temperature varies across the satellite. Expression |68| must be integrated appropriately across the surface to determine 
the net thermal thrust. 

In fa ct, the thermal thrust due to the metal is negligible because the tungsten is essentially isothermal (see equation 
(122)) and therefore radiates isotropically. The axial force is therefore supplied by a sum over CCRs: 


Faxial — 


2e g i,iRcnrR 2 

3c 


tij cos 9jTj 


(69) 
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Fig. 6. Histogram of the eclipse duration (top panels), the daily average metal temperature (middle panels), and the thermal 
drag (bottom panels) as a function of number of days since launch, for 126 days, assuming the clean-glass IR absorptivity 
otgijn. = 0.82 (left column) and the dirty-glass IR absorptivity a' g ijR = 0.6. The average thermal drag over the last 120 days is 
—0.59 pm/s 2 for clean glass and —0.36 pm/s 2 for dirty glass. 


where the factor cos 9j serves to project on the spin axis. Next, to obtain the along-track component, we project the 
axial force on the direction of motion. Then we average over one orbit, and divide by the satellite mass to obtain the 
along-track acceleration: 


along—track 


1 


J q 


F a xial{t)S ■ V sa t(t, k)dt 


(70) 


where v(t, k ) is the normalized velocity vector of the satellite. On the left column of fig. |bj we show histograms of the 
thermal drag, the metal temperature and the eclipse duration for the first 126 days using the clean-glass absorptivity 
in the IR a g ijR = 0.82. The average thermal drag over the last 120 days of the first 126-day period is found to be: 


{^along-track) days — 0.59 pm/s 


(71) 
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Fig. 7. Top row: Plot of Tw (left panel) and Tccr (right panel) (sum of the constant, first and second harmonics of the orbital 
frequency) over one period, on day 0, using the clean-glass absorptivity a g i RR = 0.82. There is no eclipse on this day, and as 
a result, the metal temperature is practically constant. The average along-track acceleration on this day is —1.0 pm/s 2 . For 
comparison, the average thermal drag over the 120-day period is —0.59 pm/s 2 . This day is not included in the average. Bottom 
row: Plot of Tw (left panel) and Tccr (right panel) (sum of the constant, first and second harmonics of the orbital frequency) 
over one period, on day 0, using the dirty-glass absorptivity a' gl IR = 0.6. The average along-track acceleration on this day is 
—0.67 pm/s 2 . For comparison, the average thermal drag over the 120-day period is —0.36 pm/s 2 . 

This is around 50 % larger than the observed value of —0.4pm/s 2 . This discrepancy can probably be explained by 
the fact that the values of all the absorptivities/emissivities are not known accurately. To illustrate this, if we lower 
the absorptivity of glass in the IR to the value a' gl IR = 0.6, a still plausible “dirty glass” value, then the average 
along-track acceleration over the 120-day period is now: 

{^along-track} 120 days — 0.36 pm/s (72) 

This time we obtain a value slightly smaller than the observed one. This suggests that there are possibly small 
unmodelled effects contributing to the observed drag —0.4/pm/s 2 . The eclipse duration, metal temperature and thermal 
drag for this “dirty glass” case are plotted on the right column in fig. [6] Also, figures [7) [8] [9] and [ T 0 | give the tungsten 
temperature and the temperatures of the various CCR rows, for days 0, 30, 60, and 90 of predictions, for both the 
clean glass a g ijR = 0.82, and the dirty glass ct' glIR = 0.60 cases. To facilitate comparison between these figures, the 
range of the plots on the left panels (for the metal temperature) is always 2 Kelvins, and the range of the plots on the 
right panels (for the CCR temperatures) is always 22 Kelvins. 


8 Conclusion 

In this paper, we have demonstrated that the anomalous along-track acceleration of LARES can indeed be explained 
by a thermal effect: the anisotropic surface temperature of the satellite results in an anisotropic radiation pressure, and 
a net force imparted on the satellite. Our computation takes into account the detailed structure of the satellite. We 
linearize the Stefan-Boltzmann law, which allows for powerful Fourier series method. This is computationally much 
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Fig. 8. Top row: Plot of TV (left panel) and Tccr (right panel) (sum of the constant, first and second harmonics of the orbital 
frequency) over one period, on day 30, using the clean-glass absorptivity u g i,iR = 0.82. The region colored gray is when the 
satellite is in the eclipse. The average along-track acceleration on this day is —0.63 pm/s 2 . For comparison, the average thermal 
drag over the 120-day period is —0.59 pm/s 2 . Bottom row: Plot of TV (left panel) and Tccr (right panel) (sum of the constant, 
first and second harmonics of the orbital frequency) over one period, on day 30, using the dirty-glass absorptivity ot' g ij R = 0.6. 
The average along-track acceleration on this day is —0.37 pm/s 2 . For comparison, the average thermal drag over the 120-day 
period is —0.36 pm/s 2 . 


less intensive than a full numerical integration along the lines of m- Our results are robust and consistent, given the 
small amplitude of temporal and spatial thermal fluctuations found. 
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A Radiative heat exchange 


Consider radiation exchange between N surfaces labelled by i (i = 1,2, • ■ ■ , N) inside an enclosure. The view factor 
Fij is defined to be the fraction of radiation leaving surface i which is intercepted by surface j. The view factor can 
be computed by evaluating a double surface integral [5Tj : 


F = _ 

n A 


COS 6i cos 6L , , , . 

- dAjdAj 


(73) 


where A is the surface element of surface i, R is the length of the straight line connecting surface elements dA with 
dAj, and 9i and 9j are the angle formed by the normals of dA and dAj with the line connecting the two elements. 
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Fig. 9. Top row: Plot of TV (left panel) and Tccr (right panel) (sum of the constant, first and second harmonics of the orbital 
frequency) over one period, on day 60, using the clean-glass absorptivity a g i t iR = 0.82. There is no eclipse on this day, and as 
a result, the metal temperature is practically constant. The average along-track acceleration on this day is —0.66 pm/s 2 . For 
comparison, the average thermal drag over the 120-day period is —0.59 pm/s 2 . Bottom row: Plot of TV (left panel) and Tccr 
(right panel) (sum of the constant, first and second harmonics of the orbital frequency) over one period, on day 60, using the 
dirty-glass absorptivity a' gi IR = 0.6. The average along-track acceleration on this day is —0.43 pm/s 2 . For comparison, the 
average thermal drag over the 120-day period is —0.36 pm/s 2 . 


Fortunately, view factors can often be deduced without having to evaluate this integral. By interchanging i and j in 
the formula above, it follows that: 

AiFij = AjFji (74) 

This is known as the reciprocity relation. Moreover, we have the summation rule 

N 

E F « = 1 (75) 

j=i 

since radiation leaving surface i must be intercepted by some surface lining the enclosure. Another observation which 
will prove useful to us is that if surface i is convex, Fa = 0 since radiation originating from i cannot be intercepted 
by i. Once all view factors between pairs of surfaces are computed, we can compute the incident irradiance on surface 
i, denoted by Hi, due to radiation by all the surfaces in the enclosure: 


~ Ai ^ FjiJjAj 

3 =1 


(76) 


where Jt is the radiosity of surface i, that is, the sum of the emitted and reflected intensities: 


Ji = e,aT 4 + (1 - ei)Hi 


(77) 
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Fig. 10. Top row: Plot of Tw (left panel) and Tccr (right panel) (sum of the constant, first and second harmonics of the 
orbital frequency) over one period, on day 90, using the clean-glass absorptivity a g ijR = 0.82. The region colored gray is when 
the satellite is in the eclipse. The average along-track acceleration on this day is —0.5 pm/s 2 . For comparison, the average 
thermal drag over the 120-day period is —0.59 pm/s 2 . Bottom row: Plot of Tw (left panel) and Tccr (right panel) (sum of 
the constant, first and second harmonics of the orbital frequency) over one period, on day 90, using the dirty-glass absorptivity 
a' g i ir = 0.6. The average along-track acceleration on this day is —0.28 pm/s 2 . For comparison, the average thermal drag over 
the first 120-day period is —0.36 pm/s 2 . 


Using the reciprocity relation (74), (76) becomes: 


N 


Hi — Fij Jj 

i=i 


Substituting into (77), we obtain a linear system of equations in the radiosities: 

N 

Ji = CiC?/ 4 + (1 — £j) Fij Jj 

i=1 


(78) 


(79) 


Once the radiosities are solved for by the standard methods of linear algebra, we can form the net heat going into 
surface i , denoted by Pi, by taking the difference between the incoming and outgoing intensities: 


P t = MH - Ji) 


(80) 


or, using ([77} again, 


Pi = 


1 - e, 


iJi-oTf) 


(81) 


Finally, we specialize to the case where the enclosure has only 2 isothermal elements, and we will assume moreover 
that element 1 is a convex surface. This scenario is useful for the purposes of a lumped-capacitance analysis of the 
CCR. The system [79] in this case becomes: 


Ji = ci oTt + (1 - ei) J 2 


(82) 





























Phuc H. Nguyen, Richard Matzner: Modelling LARES temperature distribution and thermal drag 


19 


J 2 = e 2 aT 2 4 + (1 - e 2 )[F 21 Jr + (1 - F 21 )J 2 ] (83) 

where we used the fact that Fn = 0, F 12 = 1 and equation ( [75| . Solving for the net heat exchange, we find the net 
heat flux going into element 1 to be: 

Pi = Aie eff a(T± - T 4 ) (84) 

where we defined an effective emissivity e e // : 


e eff 



ei £2 


(85) 


B Numerical values of the constants 


Semi-major axis of LARES’s orbit 

a 

7810 km 

Angular radius of Earth as viewed from LARES 

Q? e 

54.55 degrees 

Absorptivity/Emissivity of CCR in the IR (clean glass) [32] 

Otgl.IR, tgl.IR 

0.82 

Absorptivity/Emissivity of CCR in the IR (dirty glass) 

a gl.IR > e gl.IR 

0.6 

Absorptivity/Emissivity of CCR in the visible 

(XgljVist €gl,vis 

0.15 

Absorptivity/Emissivity of unoxidized tungsten in the IR at 500° C 

a W,IRi ew,IR. 

0.07 

Absorptivity/Emissivity of tungsten to sunlight 

®-W,visi €\V,vis 

0.45 

Obliquity of the ecliptic 

P 

23.5 degrees 

Mean specific heat of Suprasil glass from 0°C to 500° C [33: 

Cgl 

964 J•kg " 1 • K " 1 

Specific heat capacity of W90 Ni 6 Cu2-4 sintered tungsten [32: 

Cw 

133.9 J • kg " 1 • K " 1 

Distance from tip of CCR to bottom of cavity 

d 

5 mm 

Altitude of LARES 

h 

1540 km 

Orbital inclination 

i 

70 degrees 

Thermal conductivity of Suprasil glass at 300° C [33] 

Rgl 

1.67 W • m _1 • K _1 

Thermal conductivity of sintered tungsten THA-18N [35i 

KW 

113 W • m _1 • K _1 

Thermal diffusivity of tungsten 

A w 

4.688 x 10 " 5 m 2 • s " 1 

Total mass of LARES 

m 

387.0 kg 

Mass of a CCR 

m CCR 

0.03329 kg 

Earth IR flux per unit solid angle per unit receiver area 

Nir 

71 W • m " 2 • sr _1 

Spinning frequency of LARES at launch 

II 

o 

0.546 rad/s 

Orbital frequency of LARES 

Wo 

9.13 x 10 _ 4 rad/s 

Total solar irradiance 

F 

1366 W • nr " 2 

Radius of CCR exterior surface 

R 

0.01905 m 

Radius of Earth’s IR-emitting atmosphere 

Re 

6407 km 

Radius of the satellite 

P^sat 

0.1820 nr 

Mass density of sintered tungsten THA-18N [35] 

Pw 

18000 kg • m " 3 

Stefan-Boltzmann constant 

a 

5.670 x 10 " 8 W • m - 2 • K " 4 

Period of LARES 

T 

114.7 min 


C Analytical computation 1: the temperature of a stationary spinning sphere bathed in 
sunlight 

In this appendix, we solve for the temperature inside a stationary, spinning metal sphere bathed in sunlight. We will 
orient the z-axis along the spin axis, and suppose that sunlight arrives at an angle /3q from the spin axis. The heat 
equation reads: 

8T 

— = \ w \7 2 T (86) 

where A w is the thermal diffusivity, defined as follows: 

(87) 

CwPW 

The boundary condition reads: 

+ e WtIR aT{r = R) 4 (88) 

r—R 


r dT 

(Xvis-L — ^ "TT 

or 
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where I is the intensity of radiation incident on the metal ball. It is straightforward to show that the intensity incident 
on the surface element at P = ( 6 , 0) takes the form: 

I(t) = $ cos 7 (f) 6 >(cos 7 (f)) (89) 

where $ is the solar constant, 0 is the Heaviside step function, 7 is the (time-dependent) angle formed by the position 
vector of P with the direction of the Sun, and 

cos (7 (t)) = cos / 3q cos 9 + sin 0 0 sin 9 cos (0 — tot) (90) 

Notice that the boundary condition is nonlinear in T due to the Stefan-Boltzmann term. As in the rest of the paper, 
we linearlize the boundary condition by letting T = T 0 + AT where T 0 is the average temperature of the sphere, and 
AT is some small deviation from the average. The boundary value problem becomes: 

^ = A w V 2 AT (91) 

+ 4eaT^AT(r = R) (92) 

r—R 

Thanks to time periodicity, we can expand the radiant intensity in a Fourier series: 

OO 

I(9,<t>,t)= £ I n (0, 4>)e inut (93) 

n =—00 


, 7 ., _ dAT 

C^vis-L ^ q 


where to is the spin frequency (for simplicity, we are not considering any orbital motion in this appendix). The Fourier 
modes are determined by: 

pir / uj 

~~ inult dt (94) 

/— 7r/w 

To evaluate, we distinguish between 3 regions on the sphere. First, consider the case when 0 < 9 < | — 0 O . In this 
case, the surface element is permanently illuminated, and 


to r /u 
In( 9 ,<£) = — 

27T J-n/u 


uo<p r n/u 

/„($, 0 ) = — / (cos/3o cos 0 + sin 0 o sin 0 cos (0 — oot))e 
27r J-TV/UI 


t dt 


(95) 


In this case, only the terms with n = 0 and n = ±1 contribute: 

Io{9, (f>) = $cosf3ocos9 

I±i{d,4>) = ^ sin0 O sin 9e^ 1 ^ 

Next, consider the case § — /?o < 9 < § + 0q. In this case, the surface element is only illuminated part of the time: 


(96) 

(97) 


H 0+0 o)/w 

I n (9, 0) = —— / (cos0o cosd + sin0 o sin d cos (0 — tot))e 

J (d>—(bn)/LJ 


t dt 


'( 0-0 0 )/^ 


where 


0o(0, 0o) = arccos (— cot 0 O cot 9) 

In this case, infinitely many values of n contribute, and the first few are: 

$ 

Io{9, 0) = — (0 0 cosfio cos 9 + sin 0o sin /3 0 sin0) 

7T 


l±i( 0 , 0 ) = -e^ 
7 r 


sin 0o cos /?o cos 9 + - (0 O + cos 0 0 sin 0 O ) sin 0 O sin 9 


(98) 

(99) 

( 100 ) 

( 101 ) 


Finally, for the region § + 0o < 0 < 77 the surface element remains in the dark at all time, and I n (9, 0) = 0 for all n. 
Notice that Iq is nothing but the average intensity over a revolution. As expected, Iq is not a function of 0 but only 
of 9. 
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We also expand the temperature in a Fourier series (with the same frequency u>): 

oo 

AT(r,0,cj>,t)= ^ n (r,M)e inwt 

n =—oo 

Since the left-hand side must be real, we have ( AT n )* = AT_ n , and 

OO 

AT(r, 9, <M) = E 2Re{AT n e inu>t } 

n =0 

Substituting into the heat equation, we find that each mode AT n satisfies the Helmholtz equation: 

(V 2 + k 2 n )AT n = 0 

where 


( 102 ) 


K = - 


Xw 


(103) 


(104) 


(105) 


Unlike the usual Helmholtz equation, however, the square of the “wavenumber” k n in this case is imaginary. For the 
special case n = 0, we have ko = 0 and the Laplace equation: 


V 2 AT 0 = 0 


The general solution is: 


OO l 


AT 0 (r,9, cj>) = EE (Ai m r l + B lm r~ l ~ 1 )Y l m (9, 4>) 


(106) 


(107) 


1=0 m=—l 


Since the temperature has to be bounded at r = 0, we set Bi = 0. Also, we will start the sum over l at l = 1 since 
the l = m = 0 term is constant and can be absorbed into the average temperature Tq. For all other values of n, k n is 
nonzero and the general solution to Helmholtz equation is: 


AT n (r,6,4>) = EE (ai mn ji(k n r) + bi mn yi(k n r))Y l m {e, </>) 


(108) 


1=0 m=—l 


where ji and yi are the spherical Bessel functions and Y™ is the spherical harmonics. Requiring boundedness at r = 0 
again, we set b[ mn = 0. The general solution for AT is then: 


AT(r, 0 ,0, i) = E Ai m r l Yr(0,4>) + EE 2Re 


&Imnjl 


Vncu/X w re 3ni/4 ) <t>W 


(109) 


l,m n£N l,m L x 7 J 

We will next compute the first few coefficients of this series. By substituting the general solution into the boundary 
condition, we obtain for n = 0: 


a vis I o (0, <!>) ~ evT* = £ A lm ( K lR l ~ l + AeaT 3 R l )Yr {0 ,0) (110) 

l,m 

and for n/ 0: 

OO l 

u V i S In(0,<j >) = E E a i m n{^k n j[{k n R) +AeaT^ji(k n R))Y l m (0, </>) (111) 

1=0 m=—l 

where the derivative of the spherical Bessel functions can be expressed in terms of the spherical Bessel functions as: 

(21 + l)j'i(z) = Iji-i(z) - (l + l)ji+i(z) ( 112 ) 

It remains to expand I n (6,(f>) in the spherical coordinates and equate coefficients: 


I n (9,<t>) = ^2 I nlm Y i m (0, 4>) 
l,m 


(113) 
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where, using orthonormality of the spherical harmonics, we find: 

T n 7T 


Inim — 


/ n (M)^r*(M) sin OdOd4 


J<t >=o Je =o 

A remarkable simplification occurs if we consider the integral over <p: 


/»2-7T 

Iuim cx / e^ n+m 

Jo 


(114) 


(115) 


it is easily seen that I n im vanishes unless m = —n. Thus, it is enough to compute the coefficients I n i-n for l > n. 
In particular, among the coefficients Ai m , only those of the form Aio contribute. As a result, this series has no in¬ 
dependence and can be rewritten as a series in the Legendre polynomials. Among the coefficients a n i m , only those of 
the form a n i- n where l > n are nonzero. 

We now proceed to compute the first few I n i m . For n = 0, the first coefficient is Iqoo and it is given by: 

(Ps/ttcos 3 [3 <P rf+^ 

Iooo = --- 1—cos (3 / arccos (— cot (3 cot 0) cos 9 sin 9d9 


(p 


y+P 


sin/3 


'f ~P 


1 — cot 2 /3 cot 2 9 sin' 9d9 


With the change of variable u = — cot [3 cot 9, these integrals can be evaluated and the final result is: 

PyAr 

1 000 — —o— 


(116) 


(117) 


In particular, 7 0 oo is independent of both w and f3. This is expected since it is a time-averaged, spatially averaged 
intensity flux. The average temperature Tq will be found from the fact that Aqo = 0: 


Tn = 


/ „ T \ 1/4 
/ CXvis^ 000 \ 


= 443.6 K 


\ei R <T2yn) 

For the next coefficients we will only evaluate the integrals in for the particular case fio = 0: 

hw = fyA 
Ivo = ^ 

The time-independent piece of the temperature variation in then found to be: 

Uvis® „ , 5 a V i S <P 


AT 0 (r,9) = 


2(k + dew,IR&ToRsat) 


r cos 0 + 


32(2 Rk + 4ew,i r&Tq R 2 ) 


r 2 (3 cos 2 6 — 1) 


Or, numerically: 


AT 0 (r,9) = 0.494-—cos 0 + 0.077 
Rsat 


R S 


(3 cos 2 0 — 1) 


(118) 

(119) 

( 120 ) 

( 121 ) 

( 122 ) 


Next, we compute the coefficients ai mn of the time-dependent piece. For n = 1 (and general /3q), the lowest order term 
is: 


J n-i = R 7 \/ y sin/3 0 (8-9sin/3o-sin3/3o) 


<P 

24 

<P 

2 


i+Po 


sin/3 0 


\ Po 


arccos (— cot /3q cot 0) sin 3 9d9 


(123) 


As a consistency check, notice that if we set (3q = 0 in the expression above, we find /n_i = 0. This makes sense since, 
for f3 0 = 0, sunlight is incident from the same direction as the spin axis, and the angle of incidence for each CCR row 
is time-independent. Thus one can expect the coefficient above to vanish. 
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D Analytical computation 2: Going beyond the linearized boundary condition 


Throughout the paper, we have systematically linearized the radiative boundary condition. In this appendix, we 
develop an approximation scheme to go beyond the linearized boundary condition. To avoid complications due to 
time-dependence, we consider a non-spinning sphere bathing in sunlight. Without loss of generality, let sunlight arrive 
from the direction 0 = 0. Inside the ball, the Laplace equation for the metal is satisfied: 


V 2 T = 0 

At the surface, the following boundary condition is satisfied: 

dT 


Kw 


dr 


(Rfi) 


+ £w,irvT(R, 0) 4 = aw,vis# cos OO (^ — — d'j 


Since the temperature is independent of 4> by symmetry, we can expand in the Legendre polynomials: 

A, 


T(r,0) = J^ 


1=0 


x/W+1 


1 J Yi\0) 


Explicitly, to second order this reads: 


A\ 


A' 


T(r, 8) = A 0 Y"(8) + ^rY?(0) + ^r 2 Y 2 °(0 ) 


V3 


75 


(124) 


(125) 


(126) 


(127) 


The idea now is to expand both sides of the boundary condition p5| as a series in the Legendre polynomials, without 
linearizing the T 4 term. To do this, we will have to expand a product of spherical harmonics in terms of the spherical 
harmonics themselves. This is provided by the following formula: 


Z1+Z2 


lT(M)lr(M) = y 

l=\ll—l2 | 


(2Zi + 1)(2Z 2 + 1) 


(70|ZiZ 2 00> (/(to i + to 2 )|/i/ 2 toito 2 ) y™ 1+m2 (0,0) 


4tt(2/ + 1) 

where (Zto|ZiZ 2 TOito 2 ) are known as the Clebsch-Gordan coefficients. In particular. 


(128) 


i=\h-h\ * 4 ’ 


(129) 


The expansion for T 4 is then: 


Zl+Z2 


Z1+Z2 Z+Z' A A A A 

T 4 (r,0)= y y y y h h y l > r i ^+ i ^ 

y ( ' 4 7 r ' ) 3/2 /01" 1 

h,l2,l l 1 ,l' 2 =0l=\l 1 -l 2 \l’=\l' 1 -l' 2 \l" = \l-l'\ K ’ V 


x (Z0|ZiZ 2 00) 2 (Z'OIZ'Z^OO)^ (Z"O|ZZ'OO) 2 l7(0) 


(130) 


We will refer to the coefficient of r ra Y \ m in the expansion above as Ci mn . Suppose we work to second order, and 
compute all Ci mn with n < 2 and Z < 2. For example, let us compute Cqoo first. To do this, we have to consider all 
possible values of Zi, Z 2 , l\ , Z 2 , Z, l' and Z" such that Zi + Z 2 + l[ + Z 2 = 0 and Z" = 0. But there is only one possible 
assignment of values: l\ = Z 2 = l[ = 1' 2 = l = l' = Z" = 0. Therefore: 


Conn — 


A ° (oo|oooo ) 6 = A ° 


y ° 00 (4 tt) 3 / 2 ' uu|uuu '' / ( 4 7r)3/2 

Proceeding similarly for the other coefficients, we find the following expansion for T 4 (to second order in n and Z): 


(131) 


T 4 (r,0) = 


1 


(47r) 3 / 2 


A 4 + ^A 2 A 2 r 2 )Y 0 °(8) + 


AiAq 


8 A 2 A\ 


V / 3(47t) 3 /2 


0 rY?(0) + —j= 0 7 9 r 2 E,°(0) 


375 (4tt) 3 /2 


The right-hand side of (125|) can also be expanded as follows: 


a w , IR $cos 00(| - 0) = a w ,m$ (^Y 0 ° + J^Y? + ^Y° 


(132) 


(133) 
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From the boundary condition then, we obtain a system of 3 nonlinear equations in 3 unknowns (Ao, A\ and A 2 ). 
Solving for these coefficients, we find: 


T(r, 9) = 443.6 + 0.495-— cos 9 + 0.077 

JA'sat 


Rs 


(3 cos 2 9 — 1) 


(134) 


Comparing with equation (122) of the previous section, we find an almost identical result. Thus, the nonlinearity in 
the boundary condition really does not matter much. 
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